查看原文
其他

R统计第14篇-以TCGA数据案例演示Cox回归分析!不存在学不会

拾光2 医科研 2023-01-19

今天是各类统计方法R语言实现的第14期,我们主要介绍Cox回归分析。

Cox回归分析

之前讲到生存分析可以通过K-M生存曲线、logrank检验以及的Cox回归分析实现,其中K-M生存曲线、logrank检验只能针对单个二分类变量进行研究,想要探索连续变量或者多个因素时,则可使用Cox回归分析。

得到的结果中β为回归系数,exp(β)为hazard ratio, 简称HR,含义类似于RR。

HR = 1,没有影响。

HR > 1,风险增加,该因素为危险因素。

HR < 1,风险降低,该因素为保护因素。

数据整理

Cox回归分析全称为cox等比例风险回归模型。

此处使用的是TCGA肝癌的数据,数据经下载整理,此处仍需再次整理。

# 载入数据
lihc<-read.table("tcga_lihc.tsv",header = T,row.names = 1,quote = "",sep = '\t')
summary(lihc)
##              samples             sample_type.samples
##  TCGA-2V-A95S-01A:  1   Primary Tumor      :377     
##  TCGA-2Y-A9GS-01A:  1   Recurrent Tumor    :  3     
##  TCGA-2Y-A9GT-01A:  1   Solid Tissue Normal: 89     
##  TCGA-2Y-A9GU-01A:  1                               
##  TCGA-2Y-A9GV-01A:  1                               
##  TCGA-2Y-A9GW-01A:  1                               
##  (Other)         :463                               
##  age_at_initial_pathologic_diagnosis  tumor_stage.diagnoses
##  Min.   :16.00                       stage i     :212      
##  1st Qu.:52.00                       stage ii    :107      
##  Median :62.00                       stage iiia  : 80      
##  Mean   :60.26                       not reported: 34      
##  3rd Qu.:70.00                       stage iiib  : 12      
##  Max.   :90.00                       stage iiic  : 11      
##  NA's   :1                           (Other)     : 13      
##  neoplasm_histologic_grade       OS           OS.time      
##    :  8                    Min.   :0.000   Min.   :   1.0  
##  G1: 68                    1st Qu.:0.000   1st Qu.: 358.0  
##  G2:227                    Median :0.000   Median : 636.0  
##  G3:153                    Mean   :0.406   Mean   : 876.7  
##  G4: 13                    3rd Qu.:1.000   3rd Qu.:1214.5  
##                            Max.   :1.000   Max.   :3675.0  
##                            NA's   :6       NA's   :6       
##  adjacent_hepatic_tissue_inflammation_extent_type     PIK3CA      
##        :161                                       Min.   :0.2124  
##  Mild  :124                                       1st Qu.:0.9635  
##  None  :162                                       Median :1.2005  
##  Severe: 22                                       Mean   :1.2202  
##                                                   3rd Qu.:1.4643  
##                                                   Max.   :2.6090  
##                                                   NA's   :45      
##       AKT1            PTEN             MYC              TP53       
##  Min.   :1.461   Min.   :0.5392   Min.   :0.3709   Min.   :0.7284  
##  1st Qu.:2.708   1st Qu.:2.4475   1st Qu.:2.7515   1st Qu.:2.3810  
##  Median :2.963   Median :2.7105   Median :3.9235   Median :2.9070  
##  Mean   :2.976   Mean   :2.7179   Mean   :3.6971   Mean   :2.8948  
##  3rd Qu.:3.321   3rd Qu.:3.0350   3rd Qu.:4.7097   3rd Qu.:3.4730  
##  Max.   :4.786   Max.   :4.1610   Max.   :7.0430   Max.   :5.2620  
##  NA's   :45      NA's   :45       NA's   :45       NA's   :45

可以看到该数据包含样本名称,属于原发、复发还是癌旁,诊断年龄,分期,分级,OS状态,OS时间,癌旁炎症情况,还有几个癌基因/抑癌基因的log2(FPKM+1)

此处简单将有缺失值的样本删去,将癌旁和复发样本删去,将stage分为1、2、3、4,grade转化为1,2,3,4,癌旁炎症情况转换为1,2,3。当然也可考虑使用一些缺失值填充的方法,可以保留更多的样本。

lihc<-lihc[lihc$sample_type.samples=="Primary Tumor",]
lihc<-na.omit(lihc)
lihc<-lihc[,-c(1,2)]
lihc<-lihc[! lihc$tumor_stage.diagnoses =="not reported" &  !  lihc$neoplasm_histologic_grade == "" & 
             !  lihc$adjacent_hepatic_tissue_inflammation_extent_type == ""& !lihc$tumor_stage.diagnoses == "(Other)" ,]
summary(lihc)
##  age_at_initial_pathologic_diagnosis tumor_stage.diagnoses
##  Min.   :16.00                       stage i   :115       
##  1st Qu.:51.00                       stage ii  : 58       
##  Median :61.00                       stage iiia: 34       
##  Mean   :59.29                       stage iiib:  5       
##  3rd Qu.:69.00                       stage iiic:  4       
##  Max.   :84.00                       stage iii :  2       
##                                      (Other)   :  3       
##  neoplasm_histologic_grade       OS            OS.time      
##    :  0                    Min.   :0.0000   Min.   :   1.0  
##  G1: 24                    1st Qu.:0.0000   1st Qu.: 409.0  
##  G2:115                    Median :0.0000   Median : 662.0  
##  G3: 77                    Mean   :0.2805   Mean   : 960.1  
##  G4:  5                    3rd Qu.:1.0000   3rd Qu.:1386.0  
##                            Max.   :1.0000   Max.   :3675.0  
##                                                             
##  adjacent_hepatic_tissue_inflammation_extent_type     PIK3CA      
##        :  0                                       Min.   :0.3078  
##  Mild  : 93                                       1st Qu.:0.9649  
##  None  :112                                       Median :1.2170  
##  Severe: 16                                       Mean   :1.2237  
##                                                   3rd Qu.:1.5000  
##                                                   Max.   :2.3290  
##                                                                   
##       AKT1            PTEN            MYC              TP53       
##  Min.   :1.505   Min.   :1.095   Min.   :0.5372   Min.   :0.7284  
##  1st Qu.:2.663   1st Qu.:2.439   1st Qu.:2.6570   1st Qu.:2.4230  
##  Median :2.970   Median :2.704   Median :3.7650   Median :2.9220  
##  Mean   :2.946   Mean   :2.724   Mean   :3.5609   Mean   :2.8823  
##  3rd Qu.:3.326   3rd Qu.:3.041   3rd Qu.:4.4950   3rd Qu.:3.5330  
##  Max.   :4.786   Max.   :4.161   Max.   :7.0430   Max.   :5.2620  
#
##此处四期太少,归入三期
lihc$tumor_stage.diagnoses<- ifelse(lihc$tumor_stage.diagnoses == "stage i",1,
       ifelse(lihc$tumor_stage.diagnoses == "stage ii"|lihc$tumor_stage.diagnoses == "stage iia"|lihc$tumor_stage.diagnoses == "stage iib"|
                lihc$tumor_stage.diagnoses == "stage iic",2,3))

#
#G4太少,归入G3
lihc$neoplasm_histologic_grade<- ifelse(lihc$neoplasm_histologic_grade == "G1",1,
       ifelse(lihc$neoplasm_histologic_grade == "G2",2,3))

lihc$
adjacent_hepatic_tissue_inflammation_extent_type<-ifelse(lihc$adjacent_hepatic_tissue_inflammation_extent_type == "None",1,
       ifelse(lihc$adjacent_hepatic_tissue_inflammation_extent_type == "Mild",2,3))
summary(lihc)
##  age_at_initial_pathologic_diagnosis tumor_stage.diagnoses
##  Min.   :16.00                       Min.   :1.000        
##  1st Qu.:51.00                       1st Qu.:1.000        
##  Median :61.00                       Median :1.000        
##  Mean   :59.29                       Mean   :1.697        
##  3rd Qu.:69.00                       3rd Qu.:2.000        
##  Max.   :84.00                       Max.   :3.000        
##  neoplasm_histologic_grade       OS            OS.time      
##  Min.   :1.000             Min.   :0.0000   Min.   :   1.0  
##  1st Qu.:2.000             1st Qu.:0.0000   1st Qu.: 409.0  
##  Median :2.000             Median :0.0000   Median : 662.0  
##  Mean   :2.262             Mean   :0.2805   Mean   : 960.1  
##  3rd Qu.:3.000             3rd Qu.:1.0000   3rd Qu.:1386.0  
##  Max.   :3.000             Max.   :1.0000   Max.   :3675.0  
##  adjacent_hepatic_tissue_inflammation_extent_type     PIK3CA      
##  Min.   :1.000                                    Min.   :0.3078  
##  1st Qu.:1.000                                    1st Qu.:0.9649  
##  Median :1.000                                    Median :1.2170  
##  Mean   :1.566                                    Mean   :1.2237  
##  3rd Qu.:2.000                                    3rd Qu.:1.5000  
##  Max.   :3.000                                    Max.   :2.3290  
##       AKT1            PTEN            MYC              TP53       
##  Min.   :1.505   Min.   :1.095   Min.   :0.5372   Min.   :0.7284  
##  1st Qu.:2.663   1st Qu.:2.439   1st Qu.:2.6570   1st Qu.:2.4230  
##  Median :2.970   Median :2.704   Median :3.7650   Median :2.9220  
##  Mean   :2.946   Mean   :2.724   Mean   :3.5609   Mean   :2.8823  
##  3rd Qu.:3.326   3rd Qu.:3.041   3rd Qu.:4.4950   3rd Qu.:3.5330  
##  Max.   :4.786   Max.   :4.161   Max.   :7.0430   Max.   :5.2620
str(lihc)
## 'data.frame':    221 obs. of  11 variables:
##  $ age_at_initial_pathologic_diagnosis             : int  84 82 81 81 80 80 80 79 79 78 ...
##  $ tumor_stage.diagnoses                           : num  1 2 1 1 3 1 1 2 2 1 ...
##  $ neoplasm_histologic_grade                       : num  2 2 3 2 2 2 2 1 2 2 ...
##  $ OS                                              : int  0 1 1 0 1 1 0 0 0 1 ...
##  $ OS.time                                         : int  10 848 410 1168 1210 688 673 1241 387 1694 ...
##  $ adjacent_hepatic_tissue_inflammation_extent_type: num  1 1 1 1 1 2 3 1 1 1 ...
##  $ PIK3CA                                          : num  1.03 1.61 1.94 1.04 1.22 ...
##  $ AKT1                                            : num  4.79 3.88 3.02 2.7 3.59 ...
##  $ PTEN                                            : num  2.44 2.66 3.08 3.02 2.67 ...
##  $ MYC                                             : num  3.38 1.32 5.41 3.85 3.69 ...
##  $ TP53                                            : num  0.882 2.604 5.262 3.719 2.878 ...

整理完毕。

首先查看分期

##此处查看分期是否是肝癌的危险因素
library("survival")
## Warning: package 'survival' was built under R version 3.6.3
library("survminer")
## Warning: package 'survminer' was built under R version 3.6.3
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 3.6.3
## Loading required package: ggpubr
## Loading required package: magrittr
fit <- coxph(Surv(OS.time, OS) ~ tumor_stage.diagnoses, data = lihc)
print(fit)
## Call:
## coxph(formula = Surv(OS.time, OS) ~ tumor_stage.diagnoses, data = lihc)
#
##                         coef exp(coef) se(coef)     z        p
## tumor_stage.diagnoses 0.4904    1.6330   0.1485 3.303 0.000956
#
## Likelihood ratio test=10.66  on 1 df, p=0.001097
## n= 221, number of events= 62
res.sum <- summary(fit)
res.sum
## Call:
## coxph(formula = Surv(OS.time, OS) ~ tumor_stage.diagnoses, data = lihc)
#
##   n= 221, number of events= 62 
#
##                         coef exp(coef) se(coef)     z Pr(>|z|)    
## tumor_stage.diagnoses 0.4904    1.6330   0.1485 3.303 0.000956 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
##                       exp(coef) exp(-coef) lower .95 upper .95
## tumor_stage.diagnoses     1.633     0.6124     1.221     2.185
#
## Concordance= 0.609  (se = 0.038 )
## Likelihood ratio test= 10.66  on 1 df,   p=0.001
## Wald test            = 10.91  on 1 df,   p=0.001
## Score (logrank) test = 11.43  on 1 df,   p=7e-04

此处p<0.05,exp(coef) = 1.6330,表明分期是肝癌的危险因素。

三个p值的含义:似然比检验,Wald检验和Score (logrank)检验 ,三者几乎等价。对于足够大的样本量,可得到的类似的结果。对于较小的样本量,它们可能有所不同。对于小样本量,似然比检通常优选。

再看分级

library("survival")
library("survminer")

fit <- coxph(Surv(OS.time, OS) ~ neoplasm_histologic_grade, data = lihc)
print(fit)
## Call:
## coxph(formula = Surv(OS.time, OS) ~ neoplasm_histologic_grade, 
##     data = lihc)
#
##                             coef exp(coef) se(coef)    z      p
## neoplasm_histologic_grade 0.4762    1.6100   0.2044 2.33 0.0198
#
## Likelihood ratio test=5.69  on 1 df, p=0.01704
## n= 221, number of events= 62
res.sum <- summary(fit)
res.sum
## Call:
## coxph(formula = Surv(OS.time, OS) ~ neoplasm_histologic_grade, 
##     data = lihc)
#
##   n= 221, number of events= 62 
#
##                             coef exp(coef) se(coef)    z Pr(>|z|)  
## neoplasm_histologic_grade 0.4762    1.6100   0.2044 2.33   0.0198 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
##                           exp(coef) exp(-coef) lower .95 upper .95
## neoplasm_histologic_grade      1.61     0.6211     1.079     2.403
#
## Concordance= 0.618  (se = 0.033 )
## Likelihood ratio test= 5.69  on 1 df,   p=0.02
## Wald test            = 5.43  on 1 df,   p=0.02
## Score (logrank) test = 5.51  on 1 df,   p=0.02

此处p<0.05,exp(coef) = 1.6330,表明分级是肝癌的危险因素,

癌旁炎症情况

library("survival")
library("survminer")

fit <- coxph(Surv(OS.time, OS) ~ adjacent_hepatic_tissue_inflammation_extent_type, data = lihc)
print(fit)
## Call:
## coxph(formula = Surv(OS.time, OS) ~ adjacent_hepatic_tissue_inflammation_extent_type, 
##     data = lihc)
#
##                                                    coef exp(coef) se(coef)
## adjacent_hepatic_tissue_inflammation_extent_type 0.1032    1.1087   0.2006
##                                                      z     p
## adjacent_hepatic_tissue_inflammation_extent_type 0.514 0.607
#
## Likelihood ratio test=0.26  on 1 df, p=0.6102
## n= 221, number of events= 62
res.sum <- summary(fit)
res.sum
## Call:
## coxph(formula = Surv(OS.time, OS) ~ adjacent_hepatic_tissue_inflammation_extent_type, 
##     data = lihc)
#
##   n= 221, number of events= 62 
#
##                                                    coef exp(coef) se(coef)
## adjacent_hepatic_tissue_inflammation_extent_type 0.1032    1.1087   0.2006
##                                                      z Pr(>|z|)
## adjacent_hepatic_tissue_inflammation_extent_type 0.514    0.607
#
##                                                  exp(coef) exp(-coef) lower .95
## adjacent_hepatic_tissue_inflammation_extent_type     1.109      0.902    0.7483
##                                                  upper .95
## adjacent_hepatic_tissue_inflammation_extent_type     1.643
#
## Concordance= 0.551  (se = 0.039 )
## Likelihood ratio test= 0.26  on 1 df,   p=0.6
## Wald test            = 0.26  on 1 df,   p=0.6
## Score (logrank) test = 0.26  on 1 df,   p=0.6

此处p>0.05,表明癌旁炎症浸润对肝癌患者生存没有影响。

年龄

library("survival")
library("survminer")

fit <- coxph(Surv(OS.time, OS) ~ age_at_initial_pathologic_diagnosis , data = lihc)
print(fit)
## Call:
## coxph(formula = Surv(OS.time, OS) ~ age_at_initial_pathologic_diagnosis, 
##     data = lihc)
#
##                                        coef exp(coef) se(coef)     z     p
## age_at_initial_pathologic_diagnosis 0.02321   1.02349  0.01063 2.184 0.029
#
## Likelihood ratio test=5.19  on 1 df, p=0.02268
## n= 221, number of events= 62
res.sum <- summary(fit)
res.sum
## Call:
## coxph(formula = Surv(OS.time, OS) ~ age_at_initial_pathologic_diagnosis, 
##     data = lihc)
#
##   n= 221, number of events= 62 
#
##                                        coef exp(coef) se(coef)     z Pr(>|z|)  
## age_at_initial_pathologic_diagnosis 0.02321   1.02349  0.01063 2.184    0.029 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
##                                     exp(coef) exp(-coef) lower .95 upper .95
## age_at_initial_pathologic_diagnosis     1.023     0.9771     1.002     1.045
#
## Concordance= 0.569  (se = 0.039 )
## Likelihood ratio test= 5.19  on 1 df,   p=0.02
## Wald test            = 4.77  on 1 df,   p=0.03
## Score (logrank) test = 4.81  on 1 df,   p=0.03

此处p<0.05,exp(coef) = 1.02349,表明年龄是肝癌的危险因素,但是影响较小。

PIK3CA、AKT1、PTEN、MYC、TP53

library("survival")
library("survminer")

fit <- coxph(Surv(OS.time, OS) ~ PIK3CA , data = lihc)
print(fit)
## Call:
## coxph(formula = Surv(OS.time, OS) ~ PIK3CA, data = lihc)
#
##           coef exp(coef) se(coef)      z   p
## PIK3CA -0.1256    0.8819   0.3257 -0.386 0.7
#
## Likelihood ratio test=0.15  on 1 df, p=0.6992
## n= 221, number of events= 62
res.sum <- summary(fit)
res.sum
## Call:
## coxph(formula = Surv(OS.time, OS) ~ PIK3CA, data = lihc)
#
##   n= 221, number of events= 62 
#
##           coef exp(coef) se(coef)      z Pr(>|z|)
## PIK3CA -0.1256    0.8819   0.3257 -0.386      0.7
#
##        exp(coef) exp(-coef) lower .95 upper .95
## PIK3CA    0.8819      1.134    0.4658      1.67
#
## Concordance= 0.516  (se = 0.043 )
## Likelihood ratio test= 0.15  on 1 df,   p=0.7
## Wald test            = 0.15  on 1 df,   p=0.7
## Score (logrank) test = 0.15  on 1 df,   p=0.7

此处p>0.05,表明PIK2CA表达情况对肝癌患者生存没有影响。

library("survival")
library("survminer")

fit <- coxph(Surv(OS.time, OS) ~ AKT1 , data = lihc)
print(fit)
## Call:
## coxph(formula = Surv(OS.time, OS) ~ AKT1, data = lihc)
#
##          coef exp(coef) se(coef)     z     p
## AKT1 0.006963  1.006987 0.245491 0.028 0.977
#
## Likelihood ratio test=0  on 1 df, p=0.9774
## n= 221, number of events= 62
res.sum <- summary(fit)
res.sum
## Call:
## coxph(formula = Surv(OS.time, OS) ~ AKT1, data = lihc)
#
##   n= 221, number of events= 62 
#
##          coef exp(coef) se(coef)     z Pr(>|z|)
## AKT1 0.006963  1.006987 0.245491 0.028    0.977
#
##      exp(coef) exp(-coef) lower .95 upper .95
## AKT1     1.007     0.9931    0.6224     1.629
#
## Concordance= 0.473  (se = 0.042 )
## Likelihood ratio test= 0  on 1 df,   p=1
## Wald test            = 0  on 1 df,   p=1
## Score (logrank) test = 0  on 1 df,   p=1

此处p>0.05,表明AKT1表达情况对肝癌患者生存没有影响。

library("survival")
library("survminer")

fit <- coxph(Surv(OS.time, OS) ~ PTEN , data = lihc)
print(fit)
## Call:
## coxph(formula = Surv(OS.time, OS) ~ PTEN, data = lihc)
#
##          coef exp(coef) se(coef)      z     p
## PTEN -0.09515   0.90923  0.28585 -0.333 0.739
#
## Likelihood ratio test=0.11  on 1 df, p=0.7398
## n= 221, number of events= 62
res.sum <- summary(fit)
res.sum
## Call:
## coxph(formula = Surv(OS.time, OS) ~ PTEN, data = lihc)
#
##   n= 221, number of events= 62 
#
##          coef exp(coef) se(coef)      z Pr(>|z|)
## PTEN -0.09515   0.90923  0.28585 -0.333    0.739
#
##      exp(coef) exp(-coef) lower .95 upper .95
## PTEN    0.9092        1.1    0.5192     1.592
#
## Concordance= 0.507  (se = 0.041 )
## Likelihood ratio test= 0.11  on 1 df,   p=0.7
## Wald test            = 0.11  on 1 df,   p=0.7
## Score (logrank) test = 0.11  on 1 df,   p=0.7

此处p>0.05,表明PTEN表达情况对肝癌患者生存没有影响。

library("survival")
library("survminer")

fit <- coxph(Surv(OS.time, OS) ~ MYC , data = lihc)
print(fit)
## Call:
## coxph(formula = Surv(OS.time, OS) ~ MYC, data = lihc)
#
##        coef exp(coef) se(coef)     z     p
## MYC 0.04345   1.04440  0.09827 0.442 0.658
#
## Likelihood ratio test=0.2  on 1 df, p=0.6577
## n= 221, number of events= 62
res.sum <- summary(fit)
res.sum
## Call:
## coxph(formula = Surv(OS.time, OS) ~ MYC, data = lihc)
#
##   n= 221, number of events= 62 
#
##        coef exp(coef) se(coef)     z Pr(>|z|)
## MYC 0.04345   1.04440  0.09827 0.442    0.658
#
##     exp(coef) exp(-coef) lower .95 upper .95
## MYC     1.044     0.9575    0.8614     1.266
#
## Concordance= 0.531  (se = 0.043 )
## Likelihood ratio test= 0.2  on 1 df,   p=0.7
## Wald test            = 0.2  on 1 df,   p=0.7
## Score (logrank) test = 0.2  on 1 df,   p=0.7

此处p>0.05,表明MYC表达情况对肝癌患者生存没有影响。

library("survival")
library("survminer")

fit <- coxph(Surv(OS.time, OS) ~ TP53 , data = lihc)
print(fit)
## Call:
## coxph(formula = Surv(OS.time, OS) ~ TP53, data = lihc)
#
##         coef exp(coef) se(coef)      z     p
## TP53 -0.2363    0.7895   0.1616 -1.462 0.144
#
## Likelihood ratio test=2.1  on 1 df, p=0.1469
## n= 221, number of events= 62
res.sum <- summary(fit)
res.sum
## Call:
## coxph(formula = Surv(OS.time, OS) ~ TP53, data = lihc)
#
##   n= 221, number of events= 62 
#
##         coef exp(coef) se(coef)      z Pr(>|z|)
## TP53 -0.2363    0.7895   0.1616 -1.462    0.144
#
##      exp(coef) exp(-coef) lower .95 upper .95
## TP53    0.7895      1.267    0.5752     1.084
#
## Concordance= 0.557  (se = 0.047 )
## Likelihood ratio test= 2.1  on 1 df,   p=0.1
## Wald test            = 2.14  on 1 df,   p=0.1
## Score (logrank) test = 2.14  on 1 df,   p=0.1

此处p>0.05,表明TP53表达情况对肝癌患者生存没有影响。

看来我对于肝癌的理解还不够深啊,挑选出来的基因表达均对肝癌患者生存没有影响,当然我挑出来的都是肿瘤中经常发生突变的基因,可能突变与否对于患者生存情况影响更大吧。

年龄、分期、分级纳入回归模型

library("survival")
library("survminer")

fit <- coxph(Surv(OS.time, OS) ~ age_at_initial_pathologic_diagnosis + tumor_stage.diagnoses + neoplasm_histologic_grade, data = lihc)
print(fit) 
## Call:
## coxph(formula = Surv(OS.time, OS) ~ age_at_initial_pathologic_diagnosis + 
##     tumor_stage.diagnoses + neoplasm_histologic_grade, data = lihc)
#
##                                        coef exp(coef) se(coef)     z        p
## age_at_initial_pathologic_diagnosis 0.02699   1.02735  0.01074 2.512 0.012012
## tumor_stage.diagnoses               0.55741   1.74614  0.15315 3.640 0.000273
## neoplasm_histologic_grade           0.56870   1.76596  0.20867 2.725 0.006424
#
## Likelihood ratio test=24.72  on 3 df, p=1.767e-05
## n= 221, number of events= 62
res.sum <- summary(fit)
res.sum
## Call:
## coxph(formula = Surv(OS.time, OS) ~ age_at_initial_pathologic_diagnosis + 
##     tumor_stage.diagnoses + neoplasm_histologic_grade, data = lihc)
#
##   n= 221, number of events= 62 
#
##                                        coef exp(coef) se(coef)     z Pr(>|z|)
## age_at_initial_pathologic_diagnosis 0.02699   1.02735  0.01074 2.512 0.012012
## tumor_stage.diagnoses               0.55741   1.74614  0.15315 3.640 0.000273
## neoplasm_histologic_grade           0.56870   1.76596  0.20867 2.725 0.006424
##                                        
## age_at_initial_pathologic_diagnosis *  
## tumor_stage.diagnoses               ***
## neoplasm_histologic_grade           ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
##                                     exp(coef) exp(-coef) lower .95 upper .95
## age_at_initial_pathologic_diagnosis     1.027     0.9734     1.006     1.049
## tumor_stage.diagnoses                   1.746     0.5727     1.293     2.357
## neoplasm_histologic_grade               1.766     0.5663     1.173     2.658
#
## Concordance= 0.69  (se = 0.033 )
## Likelihood ratio test= 24.72  on 3 df,   p=2e-05
## Wald test            = 23.63  on 3 df,   p=3e-05
## Score (logrank) test = 24.6  on 3 df,   p=2e-05

可以看到三个p<0.05(似然比、Wald和Score (logrank)),表明模型是显著的。同时三个协变量均p<0.05,同时HR>1,表明三者均为独立危险因素。

好了,今天的R语言实现统计方法系列推文暂时告一段落,我们下次再见吧!小伙伴们如果有什么统计上的问题,或者如果想要学习什么方面的生物信息内容,可以在微信群或者知识星球提问,没准哪天的推文就是专门解答你的问题哦!


加入全国医学硕博交流群,与更多同行交流

一个活跃的学术交流平台,已经有300多位小伙伴跟我们一起,每天研讨氛围热烈,我自己每天都有收获,加群请备注下单位专业,质量先行。

您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存